Changing distribution and abundance of small pelagics may drive changes in predator distributions, affecting predator availability to fisheries and surveys. However, small pelagic fish are difficult to survey directly, so we developed a novel method of assessing small pelagic fish aggregate abundance via predator diet data. We used piscivore diet data collected from multiple bottom trawl surveys within a Vector Autoregressive Spatio-Temporal (VAST) model to assess trends of small pelagics on the Northeast US shelf. The goal was to develop a spatial “forage index” to inform survey and/or fishery availability in the bluefish (Pomatomus saltatrix) stock assessment. Using spring and fall surveys from 1973-2020, 20 small pelagic groups were identified as major bluefish prey using the diet data. Then, predators were grouped by diet similarity to identify 19 piscivore species with the most similar diet to bluefish in the region. Diets from all 20 piscivores were combined for the 20 prey groups at each surveyed location, and the total weight of small pelagic prey per predator stomach at each location was input into a Poisson-link delta model to estimate expected prey mass per predator stomach. Best fit models included spatial and spatio-temporal random effects, with predator mean length, number of predator species, and sea surface temperature as catchability covariates. Spring and fall prey indices were split into inshore and offshore areas to reflect changing prey availability over time in areas available to the recreational fishery and the bottom trawl survey, and also to contribute to regional ecosystem reporting.
The objective of this work was to create a “prey index” to evaluate changes in bluefish prey over time and in space using Vector Autoregressive Spatio-Temporal modeling (VAST (Thorson and Barnett, 2017; Thorson, 2019)). This approach was patterned on Ng et al. (2021), which used predator stomach data to create a biomass index for a single prey, Atlantic herring. Expected biomass of herring per stomach was estimated in VAST with 2 linear predictors, the number of herring per stomach and the average weight of herring in a stomach.
We adapted the approach of Ng et al. (2021) to get an index for “bluefish prey” in aggregate rather than a single prey species. Further, we include inshore and offshore regions by combining results across regional bottom trawl surveys surveys as was done for summer flounder biomass in Perretti and Thorson (2019). Finally, since bluefish themselves are somewhat sparsely sampled by the surveys, we aggregate all predators that have a similar diet composition to bluefish to better represent bluefish prey biomass.
We characterize mean weight of bluefish prey from all piscivores caught at each survey location and model that over time/space. Covariates potentially affecting perceived patterns in the bluefish prey index include number of predators, size composition of predators, and sea surface temperature (SST) at each survey location.
Therefore, the steps involved to estimate the forage index included defining the input dataset, and running multiple configurations of the VAST model. Steps involved in defining the dataset included defining “bluefish prey”, defining a set of piscivore predators with similar diets to bluefish, integrating diet data from two regional surveys, and integrating supplementary SST data to fill gaps in in-situ temperature data measurements. Steps involved in running the VAST model included decisions on spatial footprint, model structure, model selection to determine if spatial and spatio-temporal random effects were supported by the data, and further model selection to determine which catchability covariates were best supported by the data. Finally, subsets of the spatial domain were defined to match bluefish assessment inputs (survey and recreational fishery CPUE) for potential use as covariates in bluefish stock assessment models, and a bias-corrected (Thorson and Kristensen, 2016) forage index for each spatial subset was generated.
This approach is generalizable to any spatial subset of the full VAST spatial domain, such that forage indices for each ecoregion on the Northeast US continental shelf, or for other piscivore predators can also be generated.
Fish food habits data are collected aboard several regional fishery independent surveys in the Northeast US. The longest time series of diet has been collected by the Northeast Fisheries Science Center (NEFSC) bottom trawl survey (Smith and Link, 2010) from south of Cape Hatteras, NC to the Scotian Shelf at the US/Canada border since the early 1970s, which represents the bulk of the data for this analysis. Using similar survey protocols to the NEFSC survey, the NorthEast Area Monitoring and Assessment Program (NEAMAP) survey has collected fish diet data in inshore waters from Cape Hatteras, NC to Cape Cod, MA since 2008.
Bluefish eat small pelagics that are not well sampled by bottom trawl surveys. Bluefish themselves are not well sampled by bottom trawl surveys. Nevertheless, the diet samples collected for bluefish indicate that anchovies, herrings, squids, butterfish, scup, and small hakes are important prey. See working paper #XX, as well as this web summary page link for NEFSC survey and this presentation link for NEAMAP survey bluefish diet composition summaries.
Using all sampled bluefish stomachs included in the NEFSC food habits database 1973-2020, we created a list of all pelagic nekton bluefish prey that had at least 10 observations (Table 1). Broad categories such as empty stomach, fish unidentified, osteichthes, unidentified animal remains, and blown stomach were not included in the prey list.
# object is called `allfh`
load(url("https://github.com/Laurels1/Condition/raw/master/data/allfh.RData"))
preycount <- allfh %>%
group_by(pdcomnam, pynam) %>%
summarise(count = n()) %>%
filter(pdcomnam != "") %>%
#arrange(desc(count))
pivot_wider(names_from = pdcomnam, values_from = count)
gencomlist <- allfh %>%
select(pynam, gencom2) %>%
distinct()
blueprey <- preycount %>%
filter(BLUEFISH > 9) %>%
filter(!pynam %in% c("EMPTY", "BLOWN",
"FISH", "OSTEICHTHYES",
"ANIMAL REMAINS",
"FISH SCALES")) %>%
#filter(!str_detect(pynam, "SHRIMP|CRAB")) %>%
left_join(gencomlist) %>%
filter(!gencom2 %in% c("ARTHROPODA", "ANNELIDA",
"CNIDARIA", "UROCHORDATA",
"ECHINODERMATA", "WORMS",
"BRACHIOPODA", "COMB JELLIES",
"BRYOZOA", "SPONGES",
"MISCELLANEOUS", "OTHER")) %>%
arrange(desc(BLUEFISH))
kableExtra::kable(blueprey[,c('pynam','BLUEFISH')],
col.names = c('Prey name', 'Bluefish stomachs (n)'),
caption = "Table 1. Prey identified in bluefish stomachs, NEFSC diet database, 1973-2020.")
| Prey name | Bluefish stomachs (n) |
|---|---|
| LOLIGO SP | 423 |
| ENGRAULIDAE | 407 |
| ANCHOA MITCHILLI | 321 |
| PEPRILUS TRIACANTHUS | 307 |
| CEPHALOPODA | 262 |
| ANCHOA HEPSETUS | 171 |
| ETRUMEUS TERES | 126 |
| AMMODYTES SP | 96 |
| STENOTOMUS CHRYSOPS | 69 |
| MERLUCCIUS BILINEARIS | 53 |
| ILLEX SP | 39 |
| CLUPEA HARENGUS | 37 |
| CLUPEIDAE | 30 |
| POMATOMUS SALTATRIX | 22 |
| ENGRAULIS EURYSTOLE | 18 |
| LOLIGO PEALEII | 17 |
| SCOMBER SCOMBRUS | 14 |
| PLEURONECTIFORMES | 13 |
| CYNOSCION REGALIS | 12 |
| BREVOORTIA TYRANNUS | 10 |
The choice of predators is largely intended to balance increasing sample size for modeling bluefish prey with using predators likely to be foraging similarly to bluefish. One extreme assumption would be to include only bluefish as predators, but there are relatively few bluefish samples due to incomplete availability to bottom trawl surveys. This would miss prey available to bluefish because we have not sampled bluefish adequately. The opposite extreme assumption would be to include all stomachs that contain any of the top bluefish prey, regardless of which species ate the prey. This would include predators that do not forage similarly to bluefish and might therefore “count” prey that are not actually available to bluefish due to habitat differences. The intermediate approach which selects a group of piscivores that forage similarly to bluefish both increases sample size and screens out the most dissimilar predators to bluefish. The three possibilities above represent different piscivore groupings.
For bluefish forage index modeling, we are selecting a set of predators that have high diet similarity to bluefish. Garrison and Link (2000) evaluated similarity of predator diets on the Northeast US shelf to develop foraging guilds. The Schoener similarity index (Schoener, 1970) was applied
to assess the dietary overlap, \(D_{ij}\), between predator pairs:
\[ D_{i,j} = 1 – 0.5 (∑ |p_{i,k} – p_{j,k}|) \]
where \(p_{i,k}\) = mean proportional volume of prey type \(k\) in predator \(i\) and \(p_{i,k}\) = mean proportional volume of prey type \(k\) in predator \(j\) (Garrison and Link, 2000).
They used NEFSC bottom trawl survey data 1973-1997. We are using diets from 1985-2020 to characterize the prey index. Therefore an additional 20+ years of diet information is available to assess whether predator diet similarity has changed. Diet similarity analysis has been completed (B. Smith, pers. comm.), with a table of prey similarity available via this link on the NEFSC food habits shiny app that was used to define feeding guilds based on 50 predators. We evaluated results from several clustering algorithms (see this link) to determine how robustly different sets of predators grouped with bluefish. The working group selected the piscivore list based on the “complete” clustering algorithm (Table 2).
This piscivore dataset better captured predators that feed similarly to bluefish (e.g. striped bass), and has a higher proportion of stations with bluefish prey than a dataset based on the Garrison and Link (2000) piscivore definition. We also evaluated a piscivore definition using only the predators that always cluster with bluefish no matter what clustering algorithm is applied. However, a dataset based on that limited piscivore list excluded predators highlighted by the bluefish experts striped bass and resulted in fewer overall stations than the either of the above piscivore definitions, with a lower proportion of stations with bluefish prey.
dietoverlap <- read_csv(here("datfromshiny/tgmat.2022-02-15.csv"))
d_dietoverlap <- dist(dietoverlap)
# again directly from https://cran.r-project.org/web/packages/dendextend/vignettes/Cluster_Analysis.html
hclust_methods <- c("ward.D", "single", "complete", "average", "mcquitty",
"median", "centroid", "ward.D2")
diet_dendlist <- dendlist()
for(i in seq_along(hclust_methods)) {
hc_diet <- hclust(d_dietoverlap, method = hclust_methods[i])
diet_dendlist <- dendlist(diet_dendlist, as.dendrogram(hc_diet))
}
names(diet_dendlist) <- hclust_methods
preds <- list()
for(i in 1:8) {
dendi <- diet_dendlist[[i]]
namei <- names(diet_dendlist)[i]
labels(dendi) <- paste(as.character(names(dietoverlap[-1]))[order.dendrogram(dendi)],
"(",labels(dendi),")",
sep = "")
#preds[[namei]] <- partition_leaves(dendi)[[which_node(dendi, c("35", "36", "37"))]]
preds[[namei]] <- partition_leaves(dendi)[[which_node(dendi, c("Bluefish..S(37)", "Bluefish..M(36)", "Bluefish..L(35)"))]]
}
dendi <- diet_dendlist$complete
labels(dendi) <- paste(as.character(names(dietoverlap[-1]))[order.dendrogram(dendi)],
"(",labels(dendi),")",
sep = "")
pisccomplete <- partition_leaves(dendi)[[which_node(dendi, c("Bluefish..S(37)", "Bluefish..M(36)", "Bluefish..L(35)"))]]
pisccompletedf <- data.frame("COMNAME" = toupper(str_remove(pisccomplete, "\\..*")),
"SizeCat" = str_remove(str_extract(pisccomplete, "\\..*[:upper:]+"), "\\.."),
"feedguild" = "pisccomplete")
kableExtra::kable(pisccompletedf[,1:2],
col.names = c('Predator name', "Size category"),
caption = "Table 2. Predators clustering with Bluefish, NEFSC diet database, 1973-2020.")
| Predator name | Size category |
|---|---|
| STRIPED BASS | M |
| SEA RAVEN | L |
| STRIPED BASS | L |
| BLUEFISH | S |
| WEAKFISH | M |
| FOURSPOT FLOUNDER | L |
| NORTHERN SHORTFIN SQUID | M |
| LONGFIN SQUID | M |
| LONGFIN SQUID | S |
| NORTHERN SHORTFIN SQUID | S |
| POLLOCK | L |
| WHITE HAKE | M |
| RED HAKE | L |
| SILVER HAKE | M |
| ATLANTIC HALIBUT | M |
| POLLOCK | XL |
| WHITE HAKE | L |
| SILVER HAKE | L |
| CUSK | L |
| GOOSEFISH | XL |
| BLUEFISH | L |
| GOOSEFISH | L |
| GOOSEFISH | M |
| GOOSEFISH | S |
| SPINY DOGFISH | M |
| THORNY SKATE | XL |
| ATLANTIC HALIBUT | L |
| SEA RAVEN | M |
| SEA RAVEN | S |
| ATLANTIC COD | XL |
| SPINY DOGFISH | L |
| SPOTTED HAKE | M |
| SUMMER FLOUNDER | M |
| SUMMER FLOUNDER | L |
| BLUEFISH | M |
| BUCKLER DORY | M |
Model comparisons led us to the best model fit using mean predator length, number of predator species, and SST at a survey station as catchability covariates. The model has been partitioned into several definitions of “inshore” and “offshore” for the stock assessment inputs. First we define a partition that is the MAB and GB areas only as the GOM is not relevant to the bluefish assessment (yet). This is called MABGB, while the full model is AllEPU. Within this partition,
Survey strata definitions are built into VAST already. The results presented here pertain only to survey strata; I’ll consider the 3 mile split separately in a different model run. I can’t find a way to have VAST do all the partitions at once, because the 3 mile split cuts across survey stratum boundaries.
Definitions for strata used in the script
bluefishdiet/VASTunivariate_bfp_allsurvs_lencovSST_inoffsplit.R:
# use only MAB, GB, GOM, SS EPUs
# leave out south of Cape Hatteras at Elizabeth's suggestion
# could also leave out SS?
# CHECK if these EPUs match what we use in SOE
bfstrata <- c(3020, 3050, 3080, 3110, 3140, 3170, 3200, 3230,
3260, 3290, 3320, 3350, 3380, 3410, 3440, 3450, 3460)
bfoffshore <- c(1010, 1730, 1690, 1650, 1050, 1060, 1090, 1100, 1250, 1200, 1190, 1610)
MAB <- c(1010:1080, 1100:1120, 1600:1750, 3010:3450, 3470, 3500, 3510)
GB <- c(1090, 1130:1210, 1230, 1250, 3460, 3480, 3490, 3520:3550)
GOM <- c(1220, 1240, 1260:1290, 1360:1400, 3560:3830)
SS <- c(1300:1352, 3840:3990)
MABinshore <- c(3010:3450, 3470, 3500, 3510)
GBinshore <- c(3460, 3480, 3490, 3520:3550)
MABoffshore <- c(1010:1080, 1100:1120, 1600:1750)
GBoffshore <- c(1090, 1130:1210, 1230, 1250)
AllEPU <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MAB, GB, GOM, SS)) %>%
select(stratum_number) %>%
distinct()
MABGB <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MAB, GB)) %>%
select(stratum_number) %>%
distinct()
MABGBinshore <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MABinshore, GBinshore)) %>%
select(stratum_number) %>%
distinct()
MABGBoffshore <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MABoffshore, GBoffshore)) %>%
select(stratum_number) %>%
distinct()
bfall <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(bfstrata, bfoffshore)) %>%
select(stratum_number) %>%
distinct()
bfallnot <- northwest_atlantic_grid %>%
filter(!(stratum_number %in% bfall)) %>%
select(stratum_number) %>%
distinct()
bfin <- northwest_atlantic_grid %>%
filter(stratum_number %in% bfstrata) %>%
select(stratum_number) %>%
distinct()
bfinnot <- northwest_atlantic_grid %>%
filter(!(stratum_number %in% bfin)) %>%
select(stratum_number) %>%
distinct()
bfoff <- northwest_atlantic_grid %>%
filter(stratum_number %in% bfoffshore) %>%
select(stratum_number) %>%
distinct()
bfoffnot <- northwest_atlantic_grid %>%
filter(!(stratum_number %in% bfoff)) %>%
select(stratum_number) %>%
distinct()
MABGBalbinshore <- MABGBinshore %>%
filter(!(stratum_number %in% bfstrata)) %>%
distinct()
MABGBoffshorebigin <- MABGB %>%
filter(!(stratum_number %in% MABGBalbinshore$stratum_number)) %>%
distinct()
Attempts to rewrite these strata to a subset of only the needed ones failed! The model did not converge! I need to better understand how stratum definitions affect model estimation. I thought–apparently incorrectly–that they didn’t. This configuration, where each strata set also has a complement included, results in model convergence.
Bias correction is based on (Thorson and Kristensen, 2016).
Run this script, the same used to test catchability covariates, but
now has all three covariates (predator mean length, number of predator
species, and SSTfill) included and bias.correct=TRUE:
# VAST attempt 2 univariate model as a script
# modified from https://github.com/James-Thorson-NOAA/VAST/wiki/Index-standardization
# Load packages
library(here)
library(dplyr)
library(VAST)
#Read in data, separate spring and fall, and rename columns for VAST:
# this dataset created in SSTmethods.Rmd
bluepyagg_stn <- readRDS(here::here("fhdat/bluepyagg_stn_all_OISST.rds"))
# make SST column that uses surftemp unless missing or 0
# there are 3 surftemp 0 values in the dataset, all with oisst > 15
bluepyagg_stn <- bluepyagg_stn %>%
dplyr::mutate(sstfill = ifelse((is.na(surftemp)|surftemp==0), oisst, surftemp))
#bluepyagg_stn <- bluepyagg_pred_stn
# filter to assessment years at Tony's suggestion
# code Vessel as AL=0, HB=1, NEAMAP=2
bluepyagg_stn_fall <- bluepyagg_stn %>%
#ungroup() %>%
filter(season_ng == "FALL",
year > 1984) %>%
mutate(AreaSwept_km2 = 1, #Elizabeth's code
#declon = -declon already done before neamap merge
Vessel = as.numeric(as.factor(vessel))-1
) %>%
dplyr::select(Catch_g = meanbluepywt, #use bluepywt for individuals input in example
Year = year,
Vessel,
AreaSwept_km2,
Lat = declat,
Lon = declon,
meanpisclen,
npiscsp,
#bottemp, #this leaves out many stations for NEFSC
#surftemp, #this leaves out many stations for NEFSC
oisst,
sstfill
) %>%
na.omit() %>%
as.data.frame()
bluepyagg_stn_spring <- bluepyagg_stn %>%
filter(season_ng == "SPRING",
year > 1984) %>%
mutate(AreaSwept_km2 =1, #Elizabeth's code
#declon = -declon already done before neamap merge
Vessel = as.numeric(as.factor(vessel))-1
) %>%
dplyr::select(Catch_g = meanbluepywt,
Year = year,
Vessel,
AreaSwept_km2,
Lat = declat,
Lon = declon,
meanpisclen,
npiscsp,
#bottemp, #this leaves out many stations for NEFSC
#surftemp, #this leaves out many stations for NEFSC
oisst,
sstfill
) %>%
na.omit() %>%
as.data.frame()
# Make settings (turning off bias.correct to save time for example)
# NEFSC strata limits https://github.com/James-Thorson-NOAA/VAST/issues/302
# use only MAB, GB, GOM, SS EPUs
# leave out south of Cape Hatteras at Elizabeth's suggestion
# could also leave out SS?
# CHECK if these EPUs match what we use in SOE
bfstrata <- c(3020, 3050, 3080, 3110, 3140, 3170, 3200, 3230,
3260, 3290, 3320, 3350, 3380, 3410, 3440, 3450, 3460)
bfoffshore <- c(1010, 1730, 1690, 1650, 1050, 1060, 1090, 1100, 1250, 1200, 1190, 1610)
MAB <- c(1010:1080, 1100:1120, 1600:1750, 3010:3450, 3470, 3500, 3510)
GB <- c(1090, 1130:1210, 1230, 1250, 3460, 3480, 3490, 3520:3550)
GOM <- c(1220, 1240, 1260:1290, 1360:1400, 3560:3830)
SS <- c(1300:1352, 3840:3990)
MABinshore <- c(3010:3450, 3470, 3500, 3510)
GBinshore <- c(3460, 3480, 3490, 3520:3550)
MABoffshore <- c(1010:1080, 1100:1120, 1600:1750)
GBoffshore <- c(1090, 1130:1210, 1230, 1250)
AllEPU <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MAB, GB, GOM, SS)) %>%
select(stratum_number) %>%
distinct()
MABGB <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MAB, GB)) %>%
select(stratum_number) %>%
distinct()
MABGBinshore <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MABinshore, GBinshore)) %>%
select(stratum_number) %>%
distinct()
MABGBoffshore <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(MABoffshore, GBoffshore)) %>%
select(stratum_number) %>%
distinct()
bfall <- northwest_atlantic_grid %>%
filter(stratum_number %in% c(bfstrata, bfoffshore)) %>%
select(stratum_number) %>%
distinct()
bfallnot <- northwest_atlantic_grid %>%
filter(!(stratum_number %in% bfall)) %>%
select(stratum_number) %>%
distinct()
bfin <- northwest_atlantic_grid %>%
filter(stratum_number %in% bfstrata) %>%
select(stratum_number) %>%
distinct()
bfinnot <- northwest_atlantic_grid %>%
filter(!(stratum_number %in% bfin)) %>%
select(stratum_number) %>%
distinct()
bfoff <- northwest_atlantic_grid %>%
filter(stratum_number %in% bfoffshore) %>%
select(stratum_number) %>%
distinct()
bfoffnot <- northwest_atlantic_grid %>%
filter(!(stratum_number %in% bfoff)) %>%
select(stratum_number) %>%
distinct()
MABGBalbinshore <- MABGBinshore %>%
filter(!(stratum_number %in% bfstrata)) %>%
distinct()
MABGBoffshorebigin <- MABGB %>%
filter(!(stratum_number %in% MABGBalbinshore$stratum_number)) %>%
distinct()
# configs
FieldConfig <- c(
"Omega1" = 0, # number of spatial variation factors (0, 1, AR1)
"Epsilon1" = 0, # number of spatio-temporal factors
"Omega2" = 0,
"Epsilon2" = 0
)
# Model selection options, FieldConfig default (all IID)
# Season_knots + suffix below
# _base No vessel overdispersion or length/number covariates (ensure same dataset)
# _len Predator mean length covariate
# _no Number of predator species covariate
# _lenno Predator mean length and number of predator species covariates
# _eta10 Overdispersion (vessel effect) in first linear predictor (prey encounter)
# _eta11 Overdispersion (vessel effect) in both linear predictors (prey wt)
RhoConfig <- c(
"Beta1" = 0, # temporal structure on years (intercepts)
"Beta2" = 0,
"Epsilon1" = 0, # temporal structure on spatio-temporal variation
"Epsilon2" = 0
)
# 0 off (fixed effects)
# 1 independent
# 2 random walk
# 3 constant among years (fixed effect)
# 4 AR1
OverdispersionConfig <- c("eta1"=0, "eta2"=0)
# eta1 = vessel effects on prey encounter rate
# eta2 = vessel effects on prey weight
strata.limits <- as.list(c("AllEPU" = AllEPU,
"MABGB" = MABGB,
"MABGBinshore" = MABGBinshore,
"MABGBoffshore" = MABGBoffshore,
"bfall" = bfall,
"bfallnot" = bfallnot,
"bfin" = bfin,
"bfinnot" = bfinnot,
"bfoff" = bfoff,
"bfoffnot" = bfoffnot,
"MABGBalbinshore" = MABGBalbinshore,
"MABGBoffshorebigin" = MABGBoffshorebigin))
settings = make_settings( n_x = 500,
Region = "northwest_atlantic",
#strata.limits = list('All_areas' = 1:1e5), full area
strata.limits = strata.limits,
purpose = "index2",
bias.correct = TRUE,
#use_anisotropy = FALSE,
#fine_scale = FALSE,
#FieldConfig = FieldConfig,
#RhoConfig = RhoConfig,
OverdispersionConfig = OverdispersionConfig
)
# select dataset and set directory for output
#########################################################
# Run model fall
season <- c("fall_500_lennosst_split_biascorrect")
working_dir <- here::here(sprintf("pyindex/allagg_%s/", season))
if(!dir.exists(working_dir)) {
dir.create(working_dir)
}
fit <- fit_model(
settings = settings,
Lat_i = bluepyagg_stn_fall$Lat,
Lon_i = bluepyagg_stn_fall$Lon,
t_i = bluepyagg_stn_fall$Year,
b_i = as_units(bluepyagg_stn_fall[,'Catch_g'], 'g'),
a_i = rep(1, nrow(bluepyagg_stn_fall)),
v_i = bluepyagg_stn_fall$Vessel,
Q_ik = as.matrix(bluepyagg_stn_fall[,c("npiscsp",
"meanpisclen",
"sstfill"
)]),
#Use_REML = TRUE,
working_dir = paste0(working_dir, "/"))
saveRDS(fit, file = paste0(working_dir, "/fit.rds"))
# Plot results
plot( fit,
working_dir = paste0(working_dir, "/"))
######################################################
# Run model spring
season <- c("spring_500_lennosst_split_biascorrect")
working_dir <- here::here(sprintf("pyindex/allagg_%s/", season))
if(!dir.exists(working_dir)) {
dir.create(working_dir)
}
fit = fit_model( settings = settings,
Lat_i = bluepyagg_stn_spring[,'Lat'],
Lon_i = bluepyagg_stn_spring[,'Lon'],
t_i = bluepyagg_stn_spring[,'Year'],
b_i = as_units(bluepyagg_stn_spring[,'Catch_g'], 'g'),
a_i = rep(1, nrow(bluepyagg_stn_spring)),
v_i = bluepyagg_stn_spring$Vessel,
Q_ik = as.matrix(bluepyagg_stn_spring[,c("npiscsp",
"meanpisclen",
"sstfill"
)]),
# Use_REML = TRUE,
working_dir = paste0(working_dir, "/"))
saveRDS(fit, file = paste0(working_dir, "/fit.rds"))
# Plot results
plot( fit,
working_dir = paste0(working_dir, "/"))
Make a lookup table:
# strata.limits <- as.list(c("AllEPU" = AllEPU,
# "MABGB" = MABGB,
# "MABGBinshore" = MABGBinshore,
# "MABGBoffshore" = MABGBoffshore,
# "bfall" = bfall,
# "bfallnot" = bfallnot,
# "bfin" = bfin,
# "bfinnot" = bfinnot,
# "bfoff" = bfoff,
# "bfoffnot" = bfoffnot,
# "MABGBalbinshore" = MABGBalbinshore,
# "MABGBoffshorebigin" = MABGBoffshorebigin))
stratlook <- data.frame(Stratum = c("Stratum_1",
"Stratum_2",
"Stratum_3",
"Stratum_4",
"Stratum_5",
"Stratum_6",
"Stratum_7",
"Stratum_8",
"Stratum_9",
"Stratum_10",
"Stratum_11",
"Stratum_12"),
Region = c("AllEPU",
"MABGB",
"MABGBinshore",
"MABGBoffshore",
"bfall",
"bfallnot",
"bfin",
"bfinnot",
"bfoff",
"bfoffnot",
"MABGBalbinshore",
"MABGBoffshorebigin"))
Plot individual time series (bias corrected) with standard errors:
splitoutput <- read.csv("pyindex/allagg_fall_500_lennosst_split_biascorrect/Index.csv")
splitoutput <- splitoutput %>%
left_join(stratlook)
ggplot(splitoutput, aes(x=Time, y=Estimate, colour=Region)) +
geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
geom_point()+
geom_line()+
facet_wrap(~Region, scales = "free_y") +
guides(colour = guide_legend(ncol=2)) +
#theme(legend.position = c(1, 0),
# legend.justification = c(1, 0))
theme(legend.position="none")
in2off <- splitoutput %>%
dplyr::select(Time, Region, Estimate, Std..Error.for.Estimate) %>%
tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
dplyr::group_by(Var) %>%
tidyr::pivot_wider(names_from = Region, values_from = value) %>%
dplyr::mutate(AlbInshore = MABGBalbinshore,
BlueInshore = bfin,
BlueOffshore = bfoff,
OthOffshore = MABGB - (bfoff + bfin + MABGBalbinshore),
SumMABGB = AlbInshore + BlueInshore + BlueOffshore + OthOffshore) %>%
dplyr::select(Time, AlbInshore, BlueInshore, BlueOffshore, OthOffshore, SumMABGB, MABGB) %>%
tidyr::pivot_longer(!c(Time, Var), names_to = "Region", values_to = "value") %>%
tidyr::pivot_wider(names_from = "Var", values_from = "value")
ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
geom_point()+
geom_line()+
#facet_wrap(~Region) + #+ , scales = "free_y"
#theme(legend.position = c(1, 0),
# legend.justification = c(1, 0))
ggtitle("Fall Prey Index, Mid-Atlantic and Georges Bank")
or as proportions (here proportion of MABGB index).
MABGBprop <- in2off %>%
#dplyr::filter(Region != "AllEPU") %>%
dplyr::select(Time, Region, Estimate) %>%
tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
dplyr::mutate(AlbInshoreprop = AlbInshore/MABGB,
BlueInshoreprop = BlueInshore/MABGB,
BlueOffshoreprop = BlueOffshore/MABGB,
OthOffshoreprop = OthOffshore/MABGB) %>%
tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate") %>%
dplyr::filter(Region %in% c("AlbInshoreprop", "BlueInshoreprop", "BlueOffshoreprop",
"OthOffshoreprop"))
ggplot(MABGBprop, aes(x=Time, y=Estimate, colour = Region)) +
geom_point()+
geom_line()+
#facet_wrap(~Region) + #+ , scales = "free_y"
#theme(legend.position = c(1, 0),
# legend.justification = c(1, 0))
ggtitle("Fall Prey Index as proportion of Mid-Atlantic and Georges Bank")
Fall density maps with covariates
diagplots <- c("Data_and_knots",
"Data_by_year",
"quantile_residuals",
"quantile_residuals_on_map",
"Aniso",
"center_of_gravity")
for(p in diagplots){
cat(" \n####", as.character(p)," \n")
cat(paste0(""))
cat(" \n")
}
Plot individual time series
splitoutput <- read.csv("pyindex/allagg_spring_500_lennosst_split_biascorrect/Index.csv")
splitoutput <- splitoutput %>%
left_join(stratlook)
ggplot(splitoutput, aes(x=Time, y=Estimate, colour=Region)) +
geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
geom_point()+
geom_line()+
facet_wrap(~Region, scales = "free_y") +
guides(colour = guide_legend(ncol=2)) +
#theme(legend.position = c(1, 0),
# legend.justification = c(1, 0))
theme(legend.position="none")
or just the indices from inshore (alb), inshore bluefish, offshore bluefish, and further out
in2off <- splitoutput %>%
dplyr::select(Time, Region, Estimate, Std..Error.for.Estimate) %>%
tidyr::pivot_longer(c(Estimate, Std..Error.for.Estimate), names_to = "Var") %>%
dplyr::group_by(Var) %>%
tidyr::pivot_wider(names_from = Region, values_from = value) %>%
dplyr::mutate(AlbInshore = MABGBalbinshore,
BlueInshore = bfin,
BlueOffshore = bfoff,
OthOffshore = MABGB - (bfoff + bfin + MABGBalbinshore),
SumMABGB = AlbInshore + BlueInshore + BlueOffshore + OthOffshore) %>%
dplyr::select(Time, AlbInshore, BlueInshore, BlueOffshore, OthOffshore, SumMABGB, MABGB) %>%
tidyr::pivot_longer(!c(Time, Var), names_to = "Region", values_to = "value") %>%
tidyr::pivot_wider(names_from = "Var", values_from = "value")
ggplot(in2off, aes(x=Time, y=Estimate, colour = Region)) +
geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
geom_point()+
geom_line()+
#facet_wrap(~Region) + #+ , scales = "free_y"
#theme(legend.position = c(1, 0),
# legend.justification = c(1, 0))
ggtitle("Spring Prey Index, Mid-Atlantic and Georges Bank")
or as proportions (here proportion of MABGB index).
MABGBprop <- in2off %>%
#dplyr::filter(Region != "AllEPU") %>%
dplyr::select(Time, Region, Estimate) %>%
tidyr::pivot_wider(names_from = Region, values_from = Estimate) %>%
dplyr::mutate(AlbInshoreprop = AlbInshore/MABGB,
BlueInshoreprop = BlueInshore/MABGB,
BlueOffshoreprop = BlueOffshore/MABGB,
OthOffshoreprop = OthOffshore/MABGB) %>%
tidyr::pivot_longer(!Time, names_to = "Region", values_to = "Estimate") %>%
dplyr::filter(Region %in% c("AlbInshoreprop", "BlueInshoreprop", "BlueOffshoreprop",
"OthOffshoreprop"))
ggplot(MABGBprop, aes(x=Time, y=Estimate, colour = Region)) +
geom_point()+
geom_line()+
#facet_wrap(~Region) + #+ , scales = "free_y"
#theme(legend.position = c(1, 0),
# legend.justification = c(1, 0))
ggtitle("Spring Prey Index as proportion of Mid-Atlantic and Georges Bank")
Spring density maps with covariates
diagplots <- c("Data_and_knots",
"Data_by_year",
"quantile_residuals",
"quantile_residuals_on_map",
"Aniso",
"center_of_gravity")
for(p in diagplots){
cat(" \n####", as.character(p)," \n")
cat(paste0(""))
cat(" \n")
}
All results in the respective
pyindex/allagg_fall_500_lennosst_split_biascorrect and
allagg_spring_500_lennsst_splitbiascorrect folders.
The full results are on google drive rather than github to save space.
Still to do: